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Abstract 

Monte Carlo simulations of an atomistic solid-on-solid model are used to 
study the effect of lattice misfit on the distribution of two-dimensional is- 
lands sizes as a function of coverage G in the submonolayer aggregation 
regime of epitaxial growth. Misfit promotes the detachment of atoms from 
the perimeter of large pseudomorphic islands and thus favors their disso- 
lution into smaller islands that relieve strain more efficiently. The num- 
ber density of islands composed of s atoms exhibits scaling in the form 
N S (Q) ~ 6/(s) 2 g(s/(s)) where (s) is the average island size. Unlike the 
case of homoepitaxy, a rate equation theory based on this observation leads 
to qualitatively different behavior than observed in the simulations. 
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The morphology that obtains when atoms of one material are deposited 
onto a substrate of a dissimilar material is a central concern in current efforts 
to fabricate nanostructures in situ during growth. Equilibrium considera- 
tions make clear that small, undislocated, three-dimensional islands achieve 
significant epitaxial strain relief by lattice relaxation at the island edges || . 
But such relaxation occurs at the edges of two-dimensional (2D) heteroepi- 
taxial islands as well []2j. Thus, since there is evidence that 2D islands form a 
template from which 3D island structures evolve |3j , it seems appropriate to 
focus attention on some of the kinetic aspects of heteroepitaxy even before a 
full monolayer has been deposited. To this end, we study here the evolution 
of 2D island size distributions with a simple simulation model of epitaxial 
growth suited to the case when the difference between the deposited material 
and the substrate is completely characterized by their lattice misfit. 

The computations reported here generalize a previously successful Monte 
Carlo model of homoepitaxy Hj. There, atoms are deposited at random onto 
the (001) surface sites of a simple cubic lattice (with unit lattice constant) at 
an average rate F. No vacancies or overhangs are permitted, but any surface 
atom can hop to any nearest neighbor site at a rate D exp(— uEn /k^T), 
where D = {2kBT/h)exp(—Es/kB,T) is the single adatom migration rate, 
En is an effective pair bond energy, and n = 0,1,2,3,4 is the number of 
lateral nearest neighbors before the hop occurs. In the present work, the 
deposition and zero-strain hopping rates are fixed by the choices T = 750 K, 
F = 0.1 s _1 , £'5=1.3 eV, and £^=0.3 eV so that the dimensionless ratio 
D/F = 6 • 10 5 . All results represent an average of at least 50 realizations on 
a lattice of size 300x300. 

We suppose [|5f that the principal effect of strain is to lower the barrier 
to detachment of atoms from a strained pseudomorphic island. In a mean 
field picture, one might choose the barrier reduction to be identical for all 
atoms in an island and equal to the average strain energy per atom of that 
island. But since strain relief occurs primarily at the island edges we instead 
adopt a scheme whereby the strain-driven barrier reduction depends on the 
local coordination number More precisely, we make the replacement 
En — > En—£(^/s) for the atoms of an island composed of s atoms where S(N) 
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is the energy per atom of a one-dimensional chain of N harmonically coupled 
atoms in contact with a rigid sinusoidal potential This choice exploits 
numerical results obtained from a fully 2D version of this model [J8j which 
demonstrate that the energy density of a square island is well approximated 
by summing the energy density from two orthogonal non-interacting chains of 
atoms. Figure 1 shows a log-log plot of the strain-induced Arrhenius factor 
exp[£ (s) / kftT] with model parameters chosen to reproduce the elastic and 
cohesive properties of typical semiconductors jlj]. In practice, our numerical 
work makes use of the the approximation to this quantity indicated by the 
dashed straight line. 

The main effect of strain is to promote the dissolution of large islands by 
atom detachment ||. This is most evident from a plot of the number den- 
sity of 2D islands of size s, N S (Q), at fixed coverage for different values of 
misfit (Figure 2). Note that the monomer population is not the only benefi- 
ciary of the detachment process. Dimers and other small islands form rapidly 
but do not dissociate at an accelerated rate because they are relatively un- 
strained. The size distribution thus both narrows and shifts toward smaller 
islands sizes. More generally, we find that the average island size progres- 
sively decreases as misfit increases at fixed coverage. We note in passing that, 
compared to homoepitaxy under otherwise identical conditions, this observa- 
tion implies that relatively better layer-by-layer growth of strained material 
may be achieved so long as kinetic barriers inhibit strain-driven 3D islanding 
or misfit dislocation generation. This is so because the probability to nucleate 
next-layer islands onto the surface of existing islands decreases very sharply 
as their size decreases [TUJ . 

For the case of homoepitaxy, it is well established [JTT|, 03, [13], 031 03 that 
the island size distribution at low coverage is described by 

N,{e)~-^g(±) (1) 

( s y (s) 

where g(x) is a scaling function and (s) is the average island size. To test (0]) 
for the case of heteroepitaxy, we plot N s (s) 2 /Q versus s/ (s) and ask whether 
the simulation data at different coverages all collapse onto a single curve. 
That this is indeed the case [T5J is illustrated in Figure 3 for 5% misfit and 
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10% < 6 < 25%. Similar scaling is found for other values of misfit so that, 
quite generally, the moments of the island size distribution are given by 

M n = £ s n N s ~ e(s)"- 1 / x n g{x)dx. (2) 

s J 

To make progress, we require the coverage dependence of (s). The simu- 
lations reveal that this quantity is an increasing but not particularly simple 
function of G for all values of misfit. But in the limited coverage range noted 
above, it turns out that the power law 

(s) ~ (3) 

represents the data well. The misfit-dependence of the exponent z can be ex- 
tracted directly from (0) [square symbols in Figure 4(a)] or from the coverage 
dependence of the density of all islands combined (dashed curves in Figure 
5) since, from (0), the latter quantity takes the form N(Q) = Mq ~ @ l ~ z in 
the coverage range of interest. The triangles in Figure 4(a) show that the 
two methods yield consistent results. 

To understand the observed monotonic decrease of z toward zero it is 
convenient to return to the full N(Q) curves in Figure 5. For our choice of 
deposition conditions, N(Q) grows very rapidly independent of misfit below 
about 1.5% coverage. New island formation then slows dramatically since 
the existing islands efficiently capture newly deposited adatoms. But the 
strain energy per atom increases as the islands grow and ejection of atoms 
from perimeter sites eventually ensues. Since the ejection rate increases as 
misfit increases, the nucleation rate of new islands from this source material 
increases similarly. Comparison with the zero-strain case in Figure 5 reveals 
the efficacy of this process. 

The coverage dependence of the number density of adatoms Aq(O) is 
shown in Figure 5 as well. The origin of the relative increase in this quantity 
as a function of misfit is clear from the foregoing. More interestingly, this 
quantity is seen to exhibit a power law variation 

Aq(9) ~ 6 r (4) 

in the same coverage interval where z was defined. Figure 4(b) illustrates 
the misfit dependence of the exponent r extracted from the solid curves in 
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Figure 5 jT/J. The obvious question now arises: can a simple theory be 



constructed that predicts the values for the exponents z and r? Previous 
simulation experience with aggregation-fragmentation phenomena suggests 
that rate equation theory may be adequate for this purpose [JTB[ . To this end, 



we write an evolution equation for the number density of each island species 
in the form f|T9l 



^ = 1 - K X N\ - iVi E K S N S + l2 N 2 + E lsN s (5) 
dN, 



' s 



dG 



= N^K^N,-! - K S N S ) - ls N s + j s+1 N s+1 . (s > 1) (6) 



These mean-field equations presume that only monomers are mobile and that 
islands grow and dissociate exclusively by the attachment and detachment of 
single monomers. The rate at which adatoms attach to an island of size s is 
assumed to take the form K s = Kqs p where, e.g., the exponent p = 1/2 for the 
present case of 2D compact islands |n| . Similarly, the rate at which adatoms 



detach from an island of size s takes the form j s = "fos" for all but the very 
smallest islands. The exponent v is deduced from plots similar to Figure 1 
to be v ~ 0.027 (n) J ' 2 where (n) is the average coordination number of the 
detaching species and the lattice misfit / is expressed in percent. Although 
1 < (n) < 2, we set (n) = 1 in what follows since the precise value of 
the coefficient of f 2 anyway depends on the material parameters used in the 
model calculation of E(y/~s). 

An exact equation of motion for the moments of N s follows immediately 
from (0),® and (@): 

^r = E((s + 1)" - s")(i^AW + 1 S N S ) - l2 N 2 + K x Nl + (7) 

But we are interested in a solution at large times only so it is valid to neglect 
low order moments and the last three terms on the right hand side of (0). 
This yields the approximate expression 

~ nKoNiMn+p-i - 7i7 M n+I/ _i. (8) 
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When (0) is inserted into (§), we obtain the following self-consistent solutions 
for the desired scaling exponents 



- z/) _i < v < p 
r = z(p — v) z = \ . (9) 

v > p 

Note that z is indeterminate for the case p=v. Otherwise, the values com- 
puted from (0) are plotted as dashed curves in Figure 4 for comparison with 
the Monte Carlo values. 

It is obvious that the prediction for the monomer exponent r show the same 
trend as the simulation results while the exponent z disagrees qualitatively. 
Presumably, the non-analytic behavior of (0) is an artifact of the simplifica- 
tions required to derive (g). But are these approximations also responsible for 
the disagreement between the simulations and the rate equations regarding 
the behavior of zl To test this, we solved (|5]) and (0) numerically using rate 
parameters identical to those used in the simulations. Scaling of the assumed 
form does occur-but only at larger times (than found from the simulations) 
where coalescence should be important. Be that as it may, the qualitative be- 
havior predicted by (|9|) is confirmed albeit with the discontinuities smoothed 
out (Figure 4). 

Quantitatively, our numerical integration confirms a recent prediction by 
Blackman and Marshall |2l| that r=0 and z=l/{2 — u) for the case u=p. 
Otherwise, our prediction that r=z=0 when v > p agrees with that of Ref. 
ZUj but our result that z > 1 when v < p is not consistent with these authors 



conclusion that "gelation takes place with the formation of an infinite cluster" 
in that regime. In any event, comparison with the Monte Carlo results clearly 
impel us to conclude that the failure of the rate equations to reproduce the 
exponent z extracted from the simulations is a real effect. 

In the simulations (and in reality) an atom that detaches from an island 
generally remains close to that island and re-attaches to it with a high prob- 
ability. But in the rate equation treatment, an adatom that detaches from 
an island becomes available for capture by all islands. Moreover, since larger 
islands capture more efficiently (p > 0), the average island size will grow at 



6 



an exaggerated rate if monomers are in sufficient supply. Since increasing 
misfit precisely has the effect of generating monomers (cf. Figure 5), (0) im- 
plies that z will increase over its homoepitaxy value. When v > p, atoms 
are ejected from larger islands at a higher rate than they are captured, many 
small islands form, and z decreases precipitously. Eventually, the strain- 
induced reduction in the detachment barrier exceeds the pair bond energy 
E N and all islands disintegrate to yield adatoms as the only adsorbed species 
(*-0). 

The present results can be combined usefully with the fact that rate equa- 
tion predictions for homoepitaxy are satisfied quantitatively only for very 
large values of D/F [jn], |T5|, 0J. Both results suggest that the theory is valid 
only when the number density of monomers is sufficiently small that the error 
associated with assigning equal capture efficiency to all islands of the same 
size is negligible. This observation might help guide future research directed 
to the formulation of an improved rate equation description of the present 
problem. 
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Figure Captions 



Figure 1: Log-log plot of the island size dependence of the strain-induced 
Arrhenius factor exp[S(y/s)/kBT]. The solid symbols denote values obtained 
analytically as described in the text. The straight dashed line is the power 
law approximation to it used in the simulations and the rate equations. 

Figure 2: Distribution of island sizes for different values of misfit / after 
0.15 monolayers have been deposited. Results were obtained at T = 750 K 
with F = 0.1 s _1 , E s = 1.3 eV, and En = 0.3 eV. 

Figure 3: Typical data collapse for 5% misfit. The data shown represent 
coverages from 10% to 25%. 

Figure 4: Misfit dependence of the scaling exponents r and z. Results 
are shown from the Monte Carlo simulations at 0%, 3%, 4%, 4.7%, and 5% 
misfit for 10% < 6 < 25% (squares and triangles), an analytic analysis of the 
rate equations (dashed lines), and a numerical solution of the rate equations 
(circles). The parameters v and the misfit / are related through v ~ 0.027/ 2 . 

Figure 5: Coverage dependence of the number density of adatoms Aq(O) 
(solid curves) and the number density of all other island species combined 
N(Q) (dashed curves) for different values of misfit. 
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